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Abstract. We investigate the numerical computation of Maafi 
cusp forms for the modular group corresponding to large eigen- 
values. We present Fourier coefficients of two cusp forms whose 
eigenvalues exceed r — 40000. These eigenvalues are the largest 
that have so far been found in the case of the modular group. 
They are larger than the 130millionth eigenvalue. 



1. Introduction 

To extend the classical theory of Dirichlet series with Euler products, 
Maafi [Maa49 studied non-analytic automorphic functions, nowadays 
called Maafi waveforms. They are defined in the upper half-plane, 

H = {z = x + iy; x,y real, y > 0}, 

equipped with the hyperbolic metric 

dx 2 + dy 2 



ds z 



y 2 
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Maafi waveforms / are real analytic eigenf unctions of the hyperbolic 
Laplacian, 

(1.1) (A + A)/(z) = 0. 

The Laplacian in the hyperbolic metric reads 



A = y 2 (— + 



d 2 d 2 



"dx 2 ' dy 2 ' 

and is invariant under the group of linear fractional transformations 

z i — * ^yz — aZ \ ' ; a, b, c, d real, ad — be = 1. 
cz + a 

This group is isomorphic to the group of matrices 



7= (c d)eSL(2,R)/{±l}. 



In addition, Maafi waveforms are required to satisfy the automorphy 
condition 

(1.2) f(-rz) = f(z) V 7 GT 
relative to a cofmite discrete subgroup 

T C SL(2,R)/{±1}, 

and to satisfy the bound 

(1.3) f( z )=0(y K ) fory^oo, 

uniformly in x for some positive constant K, and similarly in the other 
cusps. Maafi waveforms which vanish in all the cusps, i.e., for which 

f(z) — > as Qz — > +oo 

and analogously at the other cusps are called Maafi cusp forms. For 
references, cf . e.g. jSeToHl IRoe66l |Hej83| ITer851 |Miy89| IVen90t Ilwa95j . 

We choose the discrete group T to be the modular group, 

T = PSL(2,Z) = SL(2,Z)/{±1}. 

It is generated by two elements 

1 l\ /0 -1 
IJ ' \1 

which are isomorphic to the translation and the inversion 

z i— > z + 1, z i — ► — , 

z 
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and has a fundamental domain which can be chosen to be 

f = r\H = {z = x + iyeH; \x\ < -, \z\>l}. 
Maafi cusp forms are square integrable over the fundamental domain 

\f(z)\ 2 dfi < oo, 
where the volume element is 

dx dy 



d\i 



y 2 



The reflection symmetry of the fundamental domain T in the line x = 
implies that the Maafi waveforms can be chosen such that they fall into 
two symmetry classes, the even functions f(x + iy) = f(—x + iy), and 
the odd ones f(x+iy) = —f(—x+iy), respectively. From the definition 
of MaaB waveforms (jl.l)) . (jl.2j) . and ([1.3)1 . it follows that they can be 
expanded into Fourier series, 

(1.4) f(z) = u (y) + ^2 a n y^ K ir (2imy) cs(2imx), 

where 



b yi +ir + hy2' ir ifr^O, 
hy^+hyHny if r = 0, 



and 



cs (x) 



2 cos(x) for the even MaaB waveforms, 
2 sin(i) for the odd ones. 



K ir (x) is the K-Bessel function (see appendix |XJ) whose order is con- 
nected with the eigenvalue A by 

A = r 2 + 1 
4 

While keeping in mind that A is the true eigenvalue, we will often call 
r to be the eigenvalue instead. 

According to the Roelcke-Selberg spectral resolution of the Lapla- 
cian [Se l56l IRoe66j , its spectrum contains both a discrete and a con- 
tinuous part. The discrete part of the spectrum is spanned by the 
constant eigenfunction / an d a countable number of Maafi cusp forms 
fiifiihi--- which we take to be ordered with increasing eigenvalues 
= A < Ai < A 2 < A 3 < . . .. The continuous part of the spectrum 
A = r 2 + i > | is spanned by the Eisenstein series E(z, f +ir) which are 
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known analytically |Maa49[ lKub73| . The functions A(^+ir)E(z, | + ir) 
are even and their Fourier coefficients are given by 

6o = A(^ + ir), 6i = A(--ir), a n = ^ |||" 

c,d£Z 
cd=n 



with 



A(s) = 7r-T(s)C(2s). 

The positive eigenvalues and their associated MaaB cusp forms are not 
known analytically. Hence one has to calculate them numerically. 

References concerning this computational work in the case of the 
modu l ar group are: jHaxTTl IASD71I IHSaTTl IC Jar78l | H^5H I?TF>??21 IHB82I 
GS84, Sta8l IWin881 



Hcj91 



ICX4S91I IHunSTl ISch.911 IHH921 |Hej92al 
|Hej92bt ISte92l IHA931 Lste94l |Hej99[ IAve03| . The firs t breakthrough 
to go beyond r = 27.284 was made by Hejhal |Hej91| who computed 
the first 123 eigenvalues and 36 more in three intervals around r « 125, 
r « 250, and r w 500, respectively. He used the truncated Fourier ex- 
pansion 

M 

(1.5) f(z) = ^ a n y*K ir {2'Kny) cs(27mx) + [[e]} 

n=l 

in the automorphy condition 

m = f{-\) 

and obtained a linear system of equations, 

M 

^a„/„(z m ) = -h{z m ), z m e J 7 , 1 < m < M - 1, 

with 

I n (z) = K- xr (2imy) cs(2imx) 



n=2 



i 



- m 



■i))3^ r (27rn3f(-i))cs(27mK(-i)), 

where, for suitable M, the error term [[e]] is of negligible size and can 
be omitted. After multiplication of this linear system of equations with 
e~5~, it was solved for successive r values on a grid. Eigenvalues were 
found by checking whether the coefficients are multiplicative, 

(1.6) ax = 1, a mp = a m a p — am, p prime, 

with the convention am = if p does not divide m. Because I n (z m ) 
gets small for large n, the linear system of equations is unstable. 
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An attempt to get around these instabilities was carried out by Stark 
[Stall], Hejhal and Arno [M93|, and Steil [Ste92l ISte94] . They used 
the eigenvalue equations 

T m /(*) = t m f(z) 

of the Hecke operators 

ad=m 
b mod d, d>0 

see Maafi |Maa49| lMaa64j . The Hecke operators T m are self-adjoint, 
commute with the Laplacian, with the symmetry of the fundamental 
domain, and amongst each other. They are multiplicative 

T m T n f(z) = T ^/W> 

d\(m,n) 
d>0 

and their eigenvalues are connected with the Fourier coefficients of the 
Maafi cusp forms by 

d m (l\t m . 

Normalizing Maafi cusp forms according to a\ = 1, the non-linear sys- 
tem of equations 

Tpf(z) = a p f(z) 
a mp = a m a P ~ af, V prime, 

allowed Steil to compute all eigenvalues up to r = 350 (4401 even and 
4776 odd eigenfunctions) and between r = 500 and 510 (395 even and 
410 odd). He was also able to compute eigenvalues around r w 4000. 

Finally, Hejhal Hej99 found a linear stable algorithm for comput- 
ing Maafi cusp forms together with their eigenvalues. His algorithm is 
based on finite Fourier transforms and implicit automorphy, and can be 
applied to holomorphic cusp forms as well as to Maafi cusp forms. Fur- 
thermore, his algorithm can also be applied to non-arithmetic groups 
and can be extended to groups whose fundamental domain has several 
cusps |SS02j . With this algorithm, Hejhal found eigenvalues around 
r w 11000. The main obstacle to go beyond was a lack of further mem- 
ory. Our goal in the present paper will be to obtain larger eigenvalues. 
We keep the main ideas, but optimize the algorithmic procedure used 
in finding the eigenvalues. Furthermore, we make careful use of the 
memory. This enables us to compute eigenvalues up to r w 40000. 
Limitations to go beyond this are due not to lack of memory, but, 
rather to CPU time. The latter (which scales with the third power of 
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r) exceeds four weeks on a 750 MHz SUN UltraSPARC-Ill processor. 
Currently, the only "larger" MaaB cusp forms available on the numer- 
ical front are those explored by Hejhal and Strombergsson HS01J m 



their recent work with waveforms of CM-type, i.e., waveforms on con- 
gruence subgroups which arise as lifts of automorphic forms on GL(1). 



2. Hejhal's algorithm 

We make use of Hejhal's algorithm Hej99 , which uses the Fourier 
expansion (jl.5j) and the automorphy condition (jl.2j) . In the present 
paper, we restrict ourselves to the modular group T = PSL(2, Z) which 
is generated by the translation z i— > z + 1 and the inversion z i— > — -. 
There do not exist small eigenvalues 0<A = r 2 + |<| for the 
modular group; see |Roe66j . Therefore, r is real and the term uo(y) in 
the Fourier expansion of MaaB cusp forms (jl.4j) vanishes. Due to the 
exponential decay of the K-Bessel function for large arguments (jA.lj) . 
and the bound 

\a n \ < d(n)ni 

for the coefficients, see |Vig83| , where d(n) counts the number of divi- 
sors of n, the absolutely convergent Fourier expansion can be truncated 
anytime we bound y from below. Given e > 0, r, and y, we determine 
the smallest M = M(e, r, y) such that the inequalities 

2wMy > r and K ir (2irMy) < e md^(K- rr (x)) 

X 

hold. Larger y allow smaller M. In all the truncated terms, i.e. within 

oo 

ll £ }}= ^ any 1 K ir (2irny) cs(2nnx), 

n=M+l 

the K-Bessel function decays exponentially in n, and already the K- 
Bessel function of the first truncated summand is smaller than e times 
most of the K-Bessel functions in the sum of (|1.5j) . Thus, the error 
[[e]] does at most marginally exceed e. The reason why [[e]] can exceed 
e somewhat is due to the possibility that the summands in (jl.5j) can 
cancel each other and that the first few coefficients a n in the truncated 
terms may occasionally be much bigger than in (|1.5|) . 

By a finite Fourier transform, the Fourier expansion (jl.5|) is solved 
for its coefficients 

(2.1) a m y^ K ir (27imy) = — f(x + iy) cs(-2vrmx) + [[e]], 
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where X is an equidistributed set of Q numbers, 

V" f _2_ _2_ ^ 2 ^ 2.1 

L 2Q ' 2Q ' ' ' ' ' 2Q ' 2Q J ' 

with 2Q > M + m. 

By automorphy we have 

/(*) = /(**), 

where z* is the T-pullback of the point z into the fundamental domain 

z* = jz, 7 g r, z* e 

Any Maafi cusp form can thus be approximated by 
(2.2) 

M 

f(x + iy) = f(x* + iy*) = ^a n y*^ ir (27rny*)cs(2W) + [[e]], 

n=l 

where is always larger or equal than the height yo of the lowest 
points in the fundamental domain T , 



y = mxn(y) = 




effectively allowing us to replace M(s, r, y) by M = M(e, r, y ). 

Choosing y smaller than y , the T-pullback of any point into the fun- 
damental domain T makes use at least once of the inversion z ^ —-, 
possibly together with the translation z i— > z + 1. This is called implicit 
automorphy, since it guarantees the invariance f(z) = /(—-), whereas 
the condition f(z) = f(z + 1) is satisfied by the Fourier expansion. 

Making use of the implicit automorphy by replacing f(x + iy) in (J2.1j) 
with the right-hand side of (J2.2|) yields 

a m y^K ir (2irmy) 

= — ^2 ^2 a ny*^K ir (2irny*) cs(2imx*) cs(-2nmx) + [[2e]] 

^ x£X n=l 

for 1 < m < M, which is the central identity in the algorithm. With 
this identity, the coefficients be determined for all m so lone as 

y < yo is chosen such that K vr {2-Rmy) does not become too small. 

Taking 1 < m < Mo and forgetting about the error [[2s]], the set of 
equations can be rewritten as 

M 

(2.3) /.Vmnjr, y)a n = 0, m>l, 

n=l 



s 
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where the matrix V = (V mn ) is given by 
V mn (r,y) = y* K ir (2irmy)8 mn 

— — — y*? Ki r (2irny*) cs(2Trnx*) cs(— 2nmx). 

Since y can always be chosen such that Ki r {2nmy) is not too small, 
the diagonal terms in the matrix V do not vanish for large m and the 
matrix is well conditioned. This makes the algorithm stable. 

We are now looking for non-trivial solutions of ()2.3|) with 1 < n < 
M that simultaneously give the eigenvalues r and the coefficients a n . 
Trivial solutions are avoided by setting a\ — 1, cf. |Miy89| . 

Since the eigenvalues r are unknown, we discretize the r axis and 
solve for each r value on this grid, 

(2.4) Yl V ™(r, y #1 )a n = -V ml (r, y* 1 ), 1 < m < M - 1, 

n=2 

where y^ 1 < yo is chosen such that Ki r (2irmy^ 1 ) is not too small for 
1 < m < M — 1. A good value to try for y^ 1 is given by 

27cM y #1 = r. 

Hejhal |Hej99| solves ()2.4j) a second time with a different y^ 2 and checks 
whether the coefficients are independent of the choice of y. 

Some words have to be said about what we mean by solving the inho- 
mogeneous system ()2.4|) . since it may happen that there is not always 
a solution unless r is an eigenvalue. By solving a linear inhomogeneous 
system of equations 

Ax = y 

we mean that we compute 

x = A~ x y 

where A~ l is determined such that A~ l A is a diagonal matrix where 
as many diagonal elements as possible are equal to one. 

3. Some improvements 

We restructure Hejhal's algorithm Hej99] in the way it finds the 
eigenvalues. Instead of solving (|2.4j) a second time, we check whether 
the coefficients a n = a* 1 obtained actually solve ()2.3|) by computing 

Mo 
n=l 
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where y* 2 = j^y* 1 is a good choice for an independent y value. Only 
if all g m vanish simultaneously can the given r be an eigenvalue and 
the computed a n 's the Fourier coefficients of a Maafi cusp form. 

The probability of finding an r value such that all g m vanish simul- 
taneously is zero because the discrete eigenvalues are of measure zero 
in the real numbers. Therefore, we make use of the intermediate value 
theorem. We let r run through a grid of discretized r values and look 
for simultaneous changes of sign in the g m . 

It is conjectured jBOGS92l IBSS92I IBol931 ISa795l IBLS96| that the 
eigenvalues of the Laplacian for even and odd cusp forms each possess 
a spacing distribution close to that of a Poisson random process. One 
therefore expects that small spacings will occur comparably often (due 
to level clustering). In order not to miss eigenvalues which lie close 
together, we have to make sure that at least one point of the r grid 
lies between any two successive eigenvalues. On the other hand, we do 
not want to waste CPU time if there are large spacings. Therefore, we 
use an adaptive algorithm which tries to predict the next best r value 
of the grid. It is based on the observation that the coefficients a n of 
two Maafi cusp forms of successive eigenvalues must differ. Assume 
that two eigenvalues lie close together and that the coefficients of the 
two Maafi cusp forms do not differ much. Numerically then both Maafi 
cusp forms would tend to be similar - which contradicts the fact that 
different Maafi cusp forms are orthogonal to each other with respect to 
the Petersson scalar product 



Maafi cusp forms corresponding to different eigenvalues are orthogonal 
because the Laplacian is an essentially self-adjoint operator. Thus, if 
successive eigenvalues lie close together, the coefficients a n must change 
fast when varying r. In contrast, if successive eigenvalues are separated 
by large spacings numerically, it turns out that often the coefficients 
change only slowly upon varying r. Defining 



our adaptive algorithm predicts the next r value of the grid such that 
the change in the coefficients is 

Mo 




a 



\<n< Mr 




(3.1) 




n=l 
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For this prediction, the last step in the r grid together with the last 
change in the coefficients is used to extrapolate linearly the choice for 
the next r value of the grid. 

However the adaptive algorithm is not rigorous. Sometimes the pre- 
diction of the next r value fails so that it is too close or too far away 
from the previous one. A small number of small steps does not bother 
us unless the step size tends to zero. But, if the step size is too large, 
such that the left-hand side of 1)3.1 Jl exceeds 0.16, we reduce the step 
size and try again with a smaller r value. 

Compared to earlier algorithms, our adaptive one tends to miss sig- 
nificantly less eigenvalues per run. 

We are searching for simultaneous sign changes in the quantities g m . 
Once we have found such in at least half of all the g m 's, we have found 
an interval [r id,r new ] which contains an eigenvalue r with high proba- 
bility. The next step is to check whether this interval really contains 
an eigenvalue, and, if so, to find this eigenvalue by some interpolation 
or bisection. 

In fact, we use a trisection which is based on a bisection together 
with a Newtonian interpolation. One first bisects the interval [r id, T n ew] 
and re-examines the sign changes. The interval with the most is then 
divided further by Newtonian interpolation, which ensures fast conver- 
gence. In the next step of the trisection, we again examine the sign 
changes and highlight that interval which contains the most simulta- 
neous sign changes in the g m 's. If there is an eigenvalue contained 
in the successive intervals of the trisection, the number of g m 's that 
simultaneously change their sign increases from step to step in the it- 
eration until the size of the interval approaches zero and the eigenvalue 
is found. In the opposite case, the number of g m 's which simultaneously 
change their sign decreases from step to step in the iteration until one 
suspects that there is no eigenvalue contained in the interval [r id, r ne w]- 

4. Results 

After some preliminary tests of our algorithm, where we computed 
some eigenvalues of the odd symmetry class around r « 10000, see 
table n an d two larger eigenvalues r = 20000.00164526 (even) and 
r = 20000.00020183 (odd), we decided to compute two Maafi cusp 
forms corresponding to eigenvalues r « 40000, one for each symmetry 
class. For the size of the error in truncating the Fourier expansion, we 
chose e = 10 -7 , which we took also for the accuracy of our K-Bessel 
function. For this cutoff, we had to take 7395 Fourier coefficients into 
account of which 938 have prime index. Finding the two eigenvalues 
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Table 1. Eigenvalues of the Maafi cusp forms around 
r « 10000 (odd symmetry). 



10000.00203541 
10000.00469659 
10000.00735313 
10000.00773954 
10000.00805085 
10000.00947268 
10000.01012235 
10000.01102222 
10000.01373844 
10000.01460515 
10000.01610617 



together with their Maafi cusp forms took four weeks of CPU time for 
each on a 750 MHz SUN UltraSPARC-Ill processor. 1.3 GB of memory 
were needed. 

Starting at r = 40000 in the upwards direction, the first even eigen- 
value was found at r = 40000.0000916; the first odd one at r = 
40000.0001644. In tables H and El we list the first few Fourier coef- 
ficients of these two forms. We checked the accuracy of our results 
with the aid of the multiplicative relations (jl.6|) . The left-hand side 
of the multiplicative relations coincides with the right-hand side up to 
a discrepancy of size 10 -3 . This means that the coefficients are only 
accurate to three digits. This is much worse than the initially intended 
accuracy e = 10~ 7 . The reason for this loss of accuracy is that we 
have computed the eigenvalues only up to an accuracy of twelve dig- 
its. Minimal deviations of the eigenvalue r lead to big changes in the 
coefficients a n . But we cannot compute the eigenvalue much more ac- 
curately without increasing the accuracy of our K-Bessel routine and 
taking more coefficients in the Fourier expansion into account. 

A different check of the accuracy of the results can be done by com- 
puting the coefficients a n a second time, with y* 1 replaced by y* 2 . But, 
with this check, one has to be careful because the coefficients may vary 
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Table 2. The first 174 Fourier coefficients of the 
MaaB cusp form corresponding to the eigenvalue r = 
40000.0000916 (even symmetry). 



fll.,,29 O30... 58 



1 




-0. 


.2674 


1 


.2094 


1 


.4654 


-0 


.1799 


-0. 


,8607 





,4629 


-0, 


,0267 


1 


.2285 


-0. 


.4445 


-0 


.2176 


-0. 


.9213 


-0 


.7499 


-0. 


,4479 


-0. 


,6494 


1. 


,3920 


-0. 


.9675 


-1. 


,3062 


1 


.4858 


0. 


,0446 





.1485 


-0. 


,7980 


-0 


.0833 


0. 


.5559 


-0 


.2480 


0, 


.1632 


-0. 


.9070 


1, 


2010 


-0 


.2211 


0, 


,0687 


-1 


.2484 


-1. 


,1888 


-0. 


.3674 


-1, 


6986 


-1 


.1703 


0, 


.9739 


-1 


.0799 


0, 


2246 


0. 


.5687 


-0. 


.4376 





.1349 


0, 


6161 





.1795 


0, 


0661 


-1 


.4044 


-0, 


,1149 





.1169 


-0, 


,4892 





.5094 


0. 


.4282 


-0 


.3000 


0, 


.1824 





.3540 


0. 


.4871 


-0 


.3473 


0. 


,1943 


0. 


.2914 


0, 


,3524 



«59...87 




38.. .116 





,5078 


-0 


.0963 


-0 


,1024 


1 


,1918 


-0 


.5538 


-1, 


.4379 


1 


.7726 





,1861 


0, 


.7256 


-0 


.6502 


0, 


.2074 


-0 


.2637 


-0 


,3046 


1 


,1780 


-0 


.0322 


-1, 


.3269 


-0 


.6666 





.1549 


-0. 


,1701 


-0 


,9765 


0, 


.2527 


-0 


,5292 


-1, 


.1143 


-0 


.1436 


1 


.2144 





.2360 





,6285 


-0 


.3695 


1 


.1076 





.0800 


1 


.6837 


-0 


.5686 


-0 


,0916 





.1612 


-0 


.5000 


0, 


.1658 


-0 


.1113 


-0 


.5918 





,0540 





.7798 


1. 


.0058 


0, 


.1639 


-1, 


.5339 





.8569 


0, 


.9039 





,2205 





.6724 


-0 


,2505 


-1 


.6999 





.9363 





,0624 


1 


,7346 


-0 


,4514 


0, 


.2351 


1. 


.4526 


-1, 


.7254 


-0. 


,0524 





,1349 



«117...145 


Ol46...174 





,2399 


1 


,3396 





,6143 





,0787 


0, 


.2756 


0, 


.6445 


0, 


.1435 


— 1, 


.0355 


-0. 


.9778 


-0. 


.1108 


-0 


.6697 


-0 


.6101 


-0 


,1000 





,7014 


0, 


.6785 


0, 


.3556 


— 0. 


.6027 


-0 


.1347 


0, 


,8777 


1 


,8005 


1 


.6805 


0, 


.0207 


1 


,1116 


-0 


.5205 


-0 


.2161 


1 


.2166 


-0 


,3686 





,0880 


1 


.5360 


— 1, 


.0575 


-0 


.0123 


1 


.0532 





,8099 


1 


.0933 


-0. 


.8063 


-0. 


.8840 


0. 


.4349 





.2573 





,2387 


-0 


.0329 


1. 


.3926 


-2, 


.0562 


0, 


.3056 


-0 


.3518 


-1, 


.0796 


-0 


.0877 


-0 


,4266 


-0 


,9383 


-0 


.1751 


-0 


.5461 


1, 


,4690 


1 


,0451 


-0 


,0367 





,5561 


1, 


.2081 


-0 


.4913 


0, 


,3580 


-0 


,0634 
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Table 3. The first 174 Fourier coefficients of the 
MaaB cusp form corresponding to the eigenvalue r = 
40000.0001644 (odd symmetry). 



Ol...29 


«30...58 


059...87 


(l; 


38.. .116 


(2117...145 


&146... 174 


1 


-0. 


.1394 





,9151 


1. 


.2840 


-0 


,1093 





,6031 


-0.5454 


0. 


.7975 


-0 
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Figure 1. Statistics of the first 938 prime Fourier co- 
efficients of the MaaB cusp form corresponding to the 
eigenvalue r = 40000.0000916. In the left figure the dis- 
tribution of the prime coefficients is rendered with points. 
The solid line is the conjectured semicircle. In the right 
figure the solid line is the integrated distribution of the 
prime coefficients, and the dashed line is the integrated 
semicircle. 



less than the size of their actual error. We did this check and found 
that the coefficients differ in the sixth digit when y^ 1 is replaced by 

y* 2 . 

All coefficients which we have computed satisfy the Ramanujan- 
Petersson conjecture 

\a p \ < 2 for all primes p. 

If the Sato- Tate conjecture is true, the prime coefficients a p of each 
MaaB cusp form are distributed according to the semicircle law 



dv(u) 



^:\/4 - u 2 du if \u\ < 2, 
otherwise. 



This means that 



lim 

N— >oo 



#{p prime; p<N} S P<N X[a,b]( a p) 



p prime 



la 



holds for any — oo < a < b < oo, where X[a,b]( u ) ^ s the indicator function 
of the interval [a, b}. The prime coefficients which we have computed 
match the Sato- Tate conjecture moderately well; see figures [T] and El 
One expects of course that the Sato- Tate conjecture is true and that the 
plots rapidly improve once one takes more coefficients (with p > M ) 
into account; cf. |HA93I ISte94] . 
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Figure 2. Statistics of the first 938 prime Fourier co- 
efficients of the MaaB cusp form corresponding to the 
eigenvalue r = 40000.0001644. In the left figure the dis- 
tribution of the prime coefficients is rendered with points. 
The solid line is the conjectured semicircle. In the right 
figure the solid line is the integrated distribution of the 
prime coefficients, and the dashed line is the integrated 
semicircle. 



5. Value distribution 

It is believed that MaaB cusp forms behave pretty much like random 
waves. In particular, in the limit of large eigenvalues, A = r 2 + | — > oo, 
a conjecture of Berry |Ber77j predicts that each MaaB cusp form has a 
Gaussian value distribution, 

1 



dp(u) 



27ra 



du, 



inside any compact regular subregion F of T . This means that 

^pj ! F X[a,b]{f{z))dLi _ 



lim 

A— >oo 



51 M 



U 



holds with variance 



1 



a' 



area(F) j F 

for any — oo<a<6<oo. Figures El and |U show the value dis- 
tribution of the MaaB cusp forms corresponding to the eigenvalues 
r = 40000.0000916 resp. r = 40000.0001644 inside a small subregion 

F = { z = x + \ y] _o.3 < x < -0.29215, 1.1 < y < 1.10785}. 

(See |Hcj99| for some analogous plots with smaller r.) 

Our numerical data agree well with Berry's conjecture, providing 
additional numerical evidence that the conjecture does hold. Plots of 




Figure 3. In the left figure the value distribution of 
the MaaB cusp form corresponding to the eigenvalue 
r = 40000.0000916 inside the region F is rendered with 
points. The solid line is the conjectured Gaussian. In the 
right figure the solid line is the integrated value distri- 
bution of the MaaB cusp form which is indistinguishable 
from the integrated Gaussian. 




Figure 4. In the left figure the value distribution of 
the MaaB cusp form corresponding to the eigenvalue 
r = 40000.0001644 inside the region F is rendered with 
points. The solid line is the conjectured Gaussian. In the 
right figure the solid line is the integrated value distri- 
bution of the MaaB cusp form which is indistinguishable 
from the integrated Gaussian. 



the two MaaB cusp forms inside the region F are given in figures El and 
Figure [7| shows the small region F inside the fundamental domain 
T. 
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Figure 5. A plot of the Maafi cusp form corresponding 
to the eigenvalue r = 40000.0000916 inside the region F. 



Appendix A. The K-Bessel function 
The K-Bessel function is defined by 



K ir (x) 



-x cosh t 



cos(rt) dt, 5fa > 0, r eC, 



see Watson |Wat44j . and is real for real arguments x and real or imag 
inary order ir. It satisfies the modified Bessel differential equation 

x 2 u"(x) + xu'(x) — (x 2 — r 2 )u(x) = 0, 
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Figure 7. The small subregion F inside the fundamen- 
tal domain T. 



MA ASS CUSP FORMS FOR LARGE EIGENVALUES 



19 



e 2 K\ 




39750 



40000 



40250 



Figure 8. K ir (x) for fixed r = 40000. 



and decays exponentially for large arguments 

71 



(A.l) 



KM ~ 




2x 



for x 



oo. 



A second linearly independent solution of the differential equation is 
the I-Bessel function 



LJx) 



' xyik 



k=0 



k\T(ir + k + 1) ' 



which grows exponentially for large arguments 



2irx 



for x — > oo. 



The amplitude of the K-Bessel function gets exponentially small if r 
increases. This can be compensated for by multiplication with the 
factor e", see figure |H1 To compute the K-Bessel function numerically 
we use asymptotic expansions for large imaginary order, r — > oo. The 
most powerful among them is the uniform asymptotic expansion 



e 2 Kir(x) ~ 22 7T 



^2 2 



oo 

E 

fc=0 



fc=0 
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where Ai(x) and Ai'(x) denote the Airy function and its derivative 
respectively. £ is defined by 

P=~, 7=-^=, -f(-0*=7" 1 -«ch- l i 9, 

and the functions v4^(x) and B^ix) are given by 

2 2k 

M s^ = (-i)'E^.(-«^^#, 

s=0 

Bkjr-jQ _ ( nW o w ^_3 a _i M 2 fc+i- 8 (7) 

r 3 s=0 

where 

An — 1, A i — — , A, = Ao_i, s ^> 2, 

' 48 48s _ ' 

and ttjfc(t) are polynomials satisfying the recursion 
u (t) = 1, 

u fc+ i(t) = V(l - t 2 )u' k (t) + lf (1 - 5t 2 H(t) dt, k>0, 

1 ° J 

see e. g. jBal66l eq. (2)], [Hal67l eqs. (18) (19) ( 20)], |(X;S91I appen- 
dix], |C4ST()2L sec tion 5], or compare with |Qlv54l eq. (4.24)], |CpU57l 
eq. (6.6)], |XS65l eas. (9.3.10) (9.3.35) (9.3.40) (9.3.41)]. All terms are 
real if x < r, and using 

2 3 

—C 2 = (-hy 1 - sec" 1 p, M-h) = (-i) fc M fc ( 7 ), 

all terms are real if x > r with positive —vy. Numerically, the uniform 
asymptotic expansion breaks down if x comes close to r. But since it 
is analytic one can expand it around the transitional point x = r and 
obtains 



e 2 K Yr (r — tr 3 ) ~ n 



/ 9 i 



A k (t) 



k=0 

oo 



1 2k 
T 3 



+ (i)W(-2*f)f>l)*^Y t small, 
r k=o rz J 



where the polynomials Ak(t) and B^it) are given in |01v52t eq. (2.42)]. 
Another useful asymptotic expansion in the transitional region is the 
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Nicholson series jMOSfifil p. 145], |Bal67l eq. (8)] 

e¥ K i {x -t x i) ^ ~ n (^) lp ( x ' *) M *)) > * small ; 

where the functions P(x,t) and Q(x,t) are defined by 

. 9 2fe 1 * \ 2 2fc i 

k=0 k=0 

and the polynomials Pk(t) and g^t) are given in |Sch54l p. 290]. Sub- 
stituting asymptotic expansions of the Airy function in the uniform 
asymptotic expansion of the K-Bessel function leads to the Hankel se- 
ries 



eir K ir {x) ~ ^/2^(sm(^{-Z)l + ^^(t) 

^ k=0 



+ cos(-(-0^ + I j^^ r « 2fe+1 (7) J, x«r 

A;=0 ' 



and the Debye series 
e^~ K; r {x) ~ 



27r(-i7)rf -exp^-^fs) ^^ k {-\i)\ x > r > 

^ fc=0 ' 

see e.g. jQEM l eqs. (2 .14) (2.19)], jSSSSI eqs. (9.7.8) (9.3.10)], [Ral67l 
eqs. (3) (5)], [CXiSHU eqs. (A9) (A10)j. Numerically, we tested all 
the given asymptotic expansions against each other to find out their 
range of applicability and their accuracy. We found by using the first 
five summands in the Hankel, the Debye and the Nicholson series, re- 
spectively, that the K-Bessel function can be approximated with an 
accuracy of at least 10 digits for r « 40000 and all x > 0. 
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